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The set of non-linear equations describing the Standard Model kinematics of the top quark an- 
tiqark production system in the dilepton decay channel has at most a four-fold ambiguity due to two 
not fully reconstructed neutrinos. Its most precise solution is of major importance for measurements 
of top quark properties like the top quark mass and tt spin correlations. Simple algebraic operations 
allow to transform the non-linear equations into a system of two polynomial equations with two 
unknowns. These two polynomials of multidegree eight can in turn be analytically reduced to one 
polynomial with one unknown by means of resultants. The obtained univariate polynomial is of 
degree sixteen. The number of its real solutions is determined analytically by means of Sturm's 
theorem, which is as well used to isolate each real solution into a unique pairwise disjoint interval. 
The solutions are polished by seeking the sign change of the polynomial in a given interval through 
binary bracketing. 



PACS numbers: PACS29.85.+C 
I. INTRODUCTION 

In 1992 R. H. Dalitz and G. R. Goldstein have published 
a numerical method based on geometrical considerations 
to solve the system of equations describingthe kinemat- 
ics of the tt decay in the dilepton channel pj . The prob- 
lem of two not fully reconstructed neutrinos - only the 
transverse components of the vector sum of their miss- 
ing energy can be measured - leads to a system of equa- 
tions which consists of as many equations as there are 
unknowns. Thus it is straight forward to solve the system 
of equations directly in contrast to a kinematic fit which 
would be appropriate in the case of an over-constrained 
problem or integration over the phase space of degrees of 
freedom in case of an under-constrained problem. Each 
of the two neutrinos contributes a two-fold ambiguity to 
the solution of the system of equations which end up to an 
over-all ambiguity of degree four. On top of those ambi- 
guities which dilute the significance of top quark property 
measurements in the dilepton channel, reconstructed ob- 
jects do typically not coincide with their corresponding 
particles which reduces the significance further. Thus it 
is not only important to solve the system of equations but 
also to compare its solutions to the particle momenta of 
simulated events. 

Next section the system of equations is introduced, fol- 
lowed by a description of the algebraic solution and its 
implementation as algorithm. Subsequently the perfor- 
mance of the numerical implementation is discussed. 



II. tt DILEPTON KINEMATICS 

The system of equations describing the kinematics of tt 
dilepton events can be expressed by the two linear and 



six non linear equations 
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The z-axis is here assumed to be parallel orientated to the 
beam axis while the x- and y-coordinates span the trans- 
verse plane. The first two equations relate the projection 
of the missing transverse energy onto one of the trans- 
verse axes (x or y) to the sum of the neutrino and an- 
tineutrino momentum components belonging to the same 
projection. The next two equations relate the energy 
of the neutrino and antineutrino, which are assumed to 
be massless in good approximation, with their momenta. 
Finally four non linear equations describe the W boson 
and top quark (antiquark) mass constraints by relating 
the invariant masses to the energy and momenta of their 
decay particles via relativistic 4- vector arithmetics. 

III. ALGEBRAIC SOLUTION 

This system of equations can be reduced to four equa- 
tions by simply substituting in the last four equations 
the neutrino and antineutrino energies by the third and 
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fourth equations and substituting the antineutrino trans- 
verse momenta by the first two equations solved to these 
momenta. In this way the four unknowns p Ux , p Vy , p Vz 
and pp z are left. One pair of equations, describing the 
t — > bW + — > frf + z^ parton branch of the event, de- 
pends on p Vz while the other pair of equations, describing 
the i — ► — > b£~D( parton branch of the event, de- 
pends on pp z . By means of ordinary algebraic operations 
both pairs can be solved to the longitudinal neutrino and 
antineutrino momentum p Vz and respectively. The 
equations can be written in the form 



p v% = a 1 ± \ I a{ + a 2 



(2) 



p Vz =h± + b 2 
for the top quark parton branch and 



p„ z = ci± \ c( + c 2 



p Pz = di ± ddj + d 



(3) 



for the anti-top quark parton branch with the coefficients 



oi = au+ ai2P„ x + aup Uv 



(4) 
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a-2 = a,2i + a 2 2P^ + a 23 p Uy + a2iP Ucc + a 2 ^p Vx p Vy + a 26 p: 

and b equivalent for the first pair of equations J2J. For 
the second pair of equations © holds analogically 

Ci = c u + c 12 p Px + c 13 p Py 

C 2 = C 2 i + C 22 pp x + C 23 Pv y + C-2AV% X + ^SPv^Pvy + C 26 pl y 

and d equivalent. The explicit expressions in terms of 
the initial equations are given in the appendix. After 
equating both equations of each pair there remain two 
equations with the two unknowns p Vx and p v . 

Again by means of ordinary algebraic operations the 
two non linear equations can be transformed into two 
polynomials of multi-degree eight. To solve these two 
polynomials to p Ux the resultant with respect to the neu- 
trino momentum p v is computed as follows. The coeffi- 
cients and monomials of the two polynomials are rewrit- 
ten in such a way that they are ordered in powers of p v 
like 



/ = hpt + hpl v + hpl + hPu v + h 



9 = 9iPt v + 92pL + 93pl + 9iPvy + 9s 



(6) 



where / and g are polynomials of the remaining un- 
knowns p Vx , p v and the coefficients f m , g n are univariate 
polynomials of p Vx . The resultant can then be obtained 



by computing the determinant of the Sylvester matrix 
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which is equated to zero. The omitted elements of the 
matrix are identical to zero. Since each element in the 
matrix is a polynomial itself the evaluation is very elab- 
orative. There are two ways to compute the determinant 
in practice. The more elegant way from a programming 
technical point of view is to invoke recursively a function 
which computes subdeterminants and consists of a very 
limited number of lines. Unfortunately it turns out that 
this approach is too time consuming. The other way is to 
let Maple compute and optimize the determinant as a 
function of the unknown p Vx and implement it. This way 
the code grows orders of magnitude in size but on the 
other hand the evaluation speeds up by orders of magni- 
tude. 

The resultant is a univariate polynomial of the form 







= hipH + kvll + h &% + h *Pll + h ^Pll + h ePlt 
+ h 7P l° + h 8 pl m + hf>p* x + h 10 pl x + hupl x (6) 
+ h 12 p 5 + h lz pl + h 14 p 3 + h 15 p 2 + h 1& p Vic + hvr 



with the remaining unknown p Vx . It is of degree 16 and 
analytical solutions of general univariate polynomials are 
only known until degree four. Abel's impossibility the- 
orem and Galois demonstrated that a univariate poly- 
nomial of degree five can in general not be solved ana- 
lytically with a finite number of additions, subtractions, 
multiplications, divisions, and root extractions [3|- Thus 
from here on the solutions of the univariate polynomial 
(|7|) have to be obtained by different means. In principle 
the problem can be reduced to an Eigenvalue problem. 
Unfortunately, in practice it turns out that the implemen- 
tation of the Eigenvalue package in Root 4j gives only 
reasonable solutions for univariate polynomials of degree 
14 and below. Finally the number of solutions is obtained 
analytically by applying Sturm's theorem [5( which con- 
sists of building a sequence of univariate polynomi- 
als h{p Vx ),b! {p Vx ),h a {p V:c ),h h {p Vx ),..., h m {p Vi ) = const., 
where hi is the first derivative of the univariate polyno- 
mial h with respect to p„ x and the following polynomials 
are the remainders of a long division of their immedi- 
ate left neighbour polynomial divided by the next left 
neighbour polynomial. The sequence ends when the last 
polynomial is a constant. In the case the constant van- 
ishes, the initial polynomial has at least one multiple real 
root which can be splitted by long division through the 
last non constant polynomial in the Sturm sequence. In 
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this case one solution is already known. The sequence 
is evaluated at two neutrino momenta p Vte (initially at 
the kinematic limits) and the difference between the num- 
ber of sign changes of the evaluated sequence at the two 
interval limits is determined. The obtained quantity cor- 
responds to the number of real solutions in the given 
interval. 

This means that the theorem of Jacques Charles 
Francois Sturm - which he has proven in 1829 Q - is 
extremely powerful since in the case of no real solutions 
no time needs to be spent for the unsuccessful attempt 
to find one. 

To reduce numerical inaccuracies, all polynomial eval- 
uations are applied using Horner's rule which factors 
out powers of the polynomial variable p Vx Further 
the solutions are separated by applying Sturm's theorem 
with varying interval boundaries. Once the solutions are 
separated in unique pairwise disjoint intervals they are 
polished by binary bracketing exploiting the knowledge 
about the sign change at the root in the given interval. 
This is possible since it is guaranteed that there is only 
one single solution in a given interval per construction 
(Now one could turn the way to solve a given Eigen- 
value problem the other way around and use the Sturm 
sequence to solve the characteristic polynomial to obtain 
the Eigenvalues). Once the solutions are found - most 
frequently there are two but never more than four (see 
fig. ^1 - they can be inserted in equations @). Such that 
these equations reduce to two univariate polynomials of 
degree four which in turn can be solved analytically to 
p v with a four fold ambiguity. The ambiguities can be 
eliminated in requiring the roots of these two polynomi- 
als to coincide since both equations have to be satisfied 
simultaneously. p Px and p Py can be simply determined 
with help of the first two equations in ffl. To determine 
the longitudinal neutrino and antineutrino momenta p Uz 
and p Pi the equations © and © can be evaluated re- 
spectively. Again the two-fold ambiguity, here due to the 
square root sign, can be resolved in requiring the solu- 
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FIG. 1: Number of solutions per event. 



tions to coincide simultaneously for both equations of one 
parton branch. 



IV. PERFORMANCE OF THE METHOD 

The univariate polynomial of p Vtt is in general very shal- 
low around zero over a broad range of neutrino momenta 
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FIG. 2: A typical univariate polynomial of degree 16 whose 
real roots in p„ x are solutions of the initial system of equations 
describing the ti dilepton kinematics. From top to bottom 
the plots are zoomed around the interesting p Vm range of the 
abscissa where two solutions are located. 
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in comparison to its maximal values in the allowed kine- 
matic range as can be concluded from the first two graphs 
in fig. El Here the kinematic range has been restricted 
to a centre of mass energy of 1.96 TeV, assuming the 
Tevatron proton anti-proton collider which has been set 
up in the Monte-Carlo event generator PYTHIA 6.220 
y| used here. Cross checks at a centre of mass energy of 
14 TeV assuming the LHC proton proton collider envi- 
ronment confirm that the performance is independent of 
particular collider settings. Only when in the graphs the 
area of the abscissa is zoomed very close to the solutions 
they can be recognised by eye. At this level the ordinate 
has already been magnified by 20 orders of magnitude. 
This explains why it is in general so difficult to find any 
solutions with numerical methods. 

In 99.9% of the events a solution can be found which is 
shown in the number of solutions per event distribution 
of fig. n] The neutrino momenta pf,° l of the solutions are 
compared to the generated ones p^ en by defining a metric 
X through 



(4) 



+ (p b c-pi: 1 ) 2 + oc-O 2 + (p^t-pi: 1 ) 2 ■ 

The solutions coincide in 99.7% of cases within real pre- 
cision to the generated neutrino momenta. Fig. [3] shows 
impressively how accurate and reliably the method is 
working. The plots in fig. 01 show the x 2 distribution 
on a linear scale. Since in practice the off-shell masses 
of the top quark and W boson resonances are not known 
the method has been applied in the following ways: The 
distribution in the first plot assumes W boson off-shell 
but top quark pole masses. It peaks at zero and its tail 
vanishes rapidly. The solution efficiency for this scenario 
amounts to 89%. The second plot assumes the pole mass 
for the top quarks and the W bosons. The number and 
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FIG. 3: Solution \ 2 denned as the difference between solved 
and generated neutrino momenta, added in quadrature, for 
the closest solution of each event. 
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FIG. 4: Minimal solution % 2 P er event. The first two plots 
differ in parton information entered into the solving proce- 
dure: a) top quark pole and W boson off-shell masses, b) top 
quark and W boson pole masses. The two lower plots show \ 2 
distributions of reconstructed events, considering both b jet 
permutations with reconstructed jets in c) and additionally 
smearing applied to jets and leptons in d). 
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Solution 


Condition 


Efficiency 


Purity 


Mean X A 


W mass known exactly 
t pole mass assumed 


0.893 


5.21 ■ 10" 5 


417.3 


t, W pole mass assumed 


0.839 





911.0 


t, W pole mass assumed 
both bb permutations 


0.890 





1166 


reconstructed 6-jets 
(parton matched) 


0.711 





3916 


wrong &-jet permutation 
(parton matched) 


0.426 





5366 


both fa-jet permutations 
(parton matched) 


0.822 





4049 


both b-jet permutations 
(parton matched, 


0.761 





4491 


jets + leptons smeared) 








both b-jet permutations 
(parton matched, 








jets + leptons smeared), 
reconstructed objects 


0.994 





2556 


100 times resolution smeared 









TABLE I: Solutions fulfilling X 2 < 10" 3 are defined as 
matched and else unmatched. The purity is determined ac- 
cording to this definition. The Mean X 2 is obtained taking 
into account matched and unmatched events. 



mean of unmatched solutions increases dramatically and 
the efficiency drops to 84%. Further an infrared-safe cone 
algorithm |9j with cone size R = 0.5 in the space spanned 
by pseudorapidity and azimuthal angle has been applied 
to the hadronic final state particles to investigate the ef- 
fect of reconstructed objects on the solutions. Requiring 
exactly two jets and two leptons and accepting the jets 
as &-tagged if they coincide within AR < 0.5 with the b 
quarks and antiquarks yields an significant degradation 
of the x 2 distribution (fig. 0]c)). The efficiency drops 
to 43% assuming the right jet quark combination. Ad- 
mitting both permutations yields an efficiency of 82%. 
The last plot has been obtained from the previous one 
in additionally smearing the leptons and jets with the 
energy resolution of the D0 detector ^(j- The \ 2 °f 
the minimal solution suffers in average another ten per- 
cent degradation and the solution efficiency drops by the 
same amount. In practice a given event passes the solving 
procedure repeatedly to improve the solution efficiency. 
Each iteration the energy of the reconstructed objects 
is randomly drawn from a probability distribution de- 
scribing the detector resolution and centred around the 
measured values. In the case of hundred such iterations 
the efficiency can be kept above 99.4% while in aver- 
age the x 2 °f the best solution decreases considerably as 
expected in comparison to solving the momenta of the 
reconstructed objects just once. 

In table [I] the efficiencies and minimal solution x 2 ' s 
are summarised. In addition the purity is given. It is 
practically zero, which means that no solutions do match 
with real precision or even merely with a x 2 better than 
10~ 3 once the off-shell masses of the top quarks and W 



bosons are not assumed to be known exactly. 

General numerical methods can compare and gauge 
their performance in terms of solution efficiency and pu- 
rity with the algebraic approach described here. 

The time consumption of the method amounts to 
about 20% of the time needed for the generation of 
the events which means if 5 • 10 6 events can be gener- 
ated in five hours an additional hour is needed to solve 
them. The strength of the method is the application of 
Sturm's theorem, such that in the case of no solutions 
the time consuming seeking and polishing of solutions 
can be saved. The bottleneck of the method is the time 
consuming evaluation of the resultant. 



V. CONCLUSIONS 

An algebraic approach to solve the ti dilepton kinemat- 
ics has been presented. The system of equations can be 
reduced to a univariate polynomial by means of resul- 
tants. The number of real roots can be determined by 
means of Sturm's theorem. Once the single roots have 
been isolated they can be polished by binary bracket- 
ing while seeking for the sign change. In this way a 
solution is found in 99.9% of cases. The solutions co- 
incide with real precision to the generated neutrino en- 
ergies and momenta in 99.7% of cases assuming that the 
reconstructible parton momenta inserted in the solving 
procedure are known exactly. Little deviations drop the 
solution efficiency considerably, at the order of tens of 
percent. In this case the solved neutrino momenta dif- 
fer already in average by the order of tens of GeV from 
the generated parton momenta. The solution efficiency 
can be re-established above the 99% level in solving a 
given event several times, varying the energy of the re- 
constructed objects each iteration randomly according 
to the energy resolution of a detector. General numeri- 
cal methods can compare their performance in terms of 
efficiency and purity to the algebraic approach whose im- 
plementation has been described here. 
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COEFFICIENTS 



Before defining the coefficients of equations and (J5J 
it is useful to introduce the following invariant masses 
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The coefficients are then given by 
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To obtain the coefficients c mn for the other parton branch 
one has to substitute W + and l + by W~ and l~ respec- 
tively. Similar holds 
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&12 
i>13 

b 2 i 
b 2 2 

&23 
b 2 4 
&25 



1 ( m t - m le+)(Pb* + Pe+) 
2{E b + E e+ ) 2 -( Pbz +p et ) 2 

fax + Pi+)0^ +Pe+) 
{E b + E e+ ) 2 - { Pbz + p e +) 2 

(Pb y +Pe+)(Pb, +P e +) 
(E b + E e+ ) 2 -( Pbz +P et ) 2 

1 4 i 4 

1 777* + 777^ + 



2m 2 m 2 M+ 



4 (E b + E e+ ) 2 - (p bz +p e +) 2 

(m 2 ~m 2 u+ ){p K +P t +) 
{E b + E e+ ) 2 -( Pbz +p et ) 2 

K-^ + )K+pJ) 

(E b + E e+ ) 2 -( Pbz +p e +) 2 
{E b + E l+ f - ( Pbv +p e +) 2 
(E b + E e+ ) 2 -( Pbz +p e +) 2 

2 (Pb* +Pe+)(Pb y +Pe+) 
(E b + E e+ ) 2 -( Pbz +p e +) 2 

(E b + E t+ ) 2 - (p by + p t +) 2 
(E b + E e+ ) 2 -( Pbz + Pe +) 2 



Again the coefficients d mn of the other parton branch 
can be obtained in substituting t, b and £ + by t, b and £~ 
respectively. The denominators are always of the type 



E 2 - p 2 z = m 2 + p\ > m 2 



Thus it is ensured that they never vanish. Running over 
1 million Monte-Carlo events does not lead to a division 
by zero. In addition, detected objects in collider exper- 
iments have always a considerable amount of transverse 
momentum which pushes the kinematics of the equations 
further away from such singularities. Therefore the theo- 
retically possible multiplication of all equations with the 
least common multiple of all denominators does not need 
to be applied. 
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